#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "mpi.h"
//mpicc -std=c99 -o cmake-build/mpi/multiplication mpi/multiplication.c && mpirun -np 4 ./cmake-build/mpi/multiplication

#define matrixARow (3)
#define matrixACol (4)
#define matrixBRow matrixACol   //A的列和B的行相等
#define matrixBCol (2)

#define matrixCRow matrixARow
#define matrixCCol matrixBCol

double matrixA[matrixARow][matrixACol];
double matrixB[matrixBRow][matrixBCol];
double matrixC[matrixCRow][matrixCCol];

int main(int argc,char*argv[])
{
    int rank,size;
    MPI_Init(&argc,&argv);
    MPI_Comm_size(MPI_COMM_WORLD,&size);
    MPI_Comm_rank(MPI_COMM_WORLD,&rank);
    
    //进程0初始化输入矩阵和结果矩阵
    if(rank==0)
    {
        matrixA[0][0]=1;
        matrixA[0][1]=2;
        matrixA[0][2]=3;
        matrixA[0][3]=4;
        matrixA[1][0]=2;
        matrixA[1][1]=3;
        matrixA[1][2]=4;
        matrixA[1][3]=5;
        matrixA[2][0]=3;
        matrixA[2][1]=4;
        matrixA[2][2]=5;
        matrixA[2][3]=6;
        
        matrixB[0][0]=1;
        matrixB[1][0]=2;
        matrixB[2][0]=3;
        matrixB[3][0]=4;
    
        matrixB[0][1]=1;
        matrixB[1][1]=2;
        matrixB[2][1]=3;
        matrixB[3][1]=4;
        
        for(int i=0;i<matrixCRow;++i)
            for(int j=0;j<matrixCCol;++j)
                matrixC[i][j]=0;
        
        for(int i=0;i<matrixARow;++i,putchar('\n'))
            for(int j=0;j<matrixACol;++j)
                printf("%6.2lf ",matrixA[i][j]);
        
        for(int i=0;i<matrixBRow;++i,putchar('\n'))
            for(int j=0;j<matrixBCol;++j)
                printf("%6.2lf ",matrixB[i][j]);
    }
    //向所有进程广播
    MPI_Bcast(matrixA[0],matrixARow*matrixACol,MPI_DOUBLE,0,MPI_COMM_WORLD);
    MPI_Bcast(matrixB[0],matrixBRow*matrixBCol,MPI_DOUBLE,0,MPI_COMM_WORLD);
    MPI_Bcast(matrixC[0],matrixCRow*matrixCCol,MPI_DOUBLE,0,MPI_COMM_WORLD);
    
    //计算
    for(int l=rank;l<matrixCRow*matrixCCol;l+=size)
    {
        int i=l/matrixCCol;
        int j=l-i*matrixCCol;
        for(int k=0;k<matrixACol;++k)
            matrixC[i][j]+=matrixA[i][k]*matrixB[k][j];
        printf("Thread %d process row %d col %d %.2lf\n",rank,i,j,matrixC[i][j]);
        fflush(stdout);
    }
    //发回结果
    MPI_Barrier(MPI_COMM_WORLD);
    MPI_Allreduce(MPI_IN_PLACE,matrixC[0],matrixCRow*matrixCCol,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
    
    //输出结果
    fflush(stdout);
    if(rank==0)
    {
        for(int i=0;i<matrixCRow;++i,putchar('\n'))
            for(int j=0;j<matrixCCol;++j)
                printf("%6.2lf ",matrixC[i][j]);
    }
    
    MPI_Finalize();
    return 0;
}
